!=================================================================================================================
 module bl_mynn
 use ccpp_kind_types,only: kind_phys

 use bl_mynn_common,only: &
        cp      , cpv       , cliq      , cice , ep_1 , ep_2  , ep_3  , grav   , karman , p1000mb , &
        r_d     , r_v       , svp1      , svp2 , svp3 , svpt0 , xlf   , xls    , xlv    , p608    , &
        t0c     , tref      , tkmin     , tv0  , gtr  , xlvcp , xlscp , rvovrd , rcp    , cphh_st , &
        cphm_st , cphh_unst , cphm_unst , b1   , b2   , zero
 use bl_mynn_subroutines


 implicit none
 private
 public:: bl_mynn_init,     &
          bl_mynn_finalize, &
          bl_mynn_run


 contains


!=================================================================================================================
 subroutine bl_mynn_init(con_cp,con_cpv,con_cice,con_cliq,con_ep1,con_ep2,con_grav,con_karman,con_p0, &
                         con_rd,con_rv,con_svp1,con_svp2,con_svp3,con_svpt0,con_xlf,con_xls,con_xlv,  &
                         errmsg,errflg)
!=================================================================================================================

!-- input arguments:
 real(kind=kind_phys),intent(in):: &
    con_cp,     &
    con_cpv,    &
    con_cice,   &
    con_cliq

 real(kind=kind_phys),intent(in):: &
    con_ep1,    &
    con_ep2

 real(kind=kind_phys),intent(in):: &
    con_grav

 real(kind=kind_phys),intent(in):: &
    con_karman

 real(kind=kind_phys),intent(in):: &
    con_p0

 real(kind=kind_phys),intent(in):: &
    con_rd,     &
    con_rv

 real(kind=kind_phys),intent(in):: &
    con_svp1,   &
    con_svp2,   &
    con_svp3,   &
    con_svpt0

 real(kind=kind_phys),intent(in):: &
    con_xlf,    &
    con_xls,    &
    con_xlv

 character(len=*),intent(out):: errmsg
 integer,intent(out):: errflg

!-----------------------------------------------------------------------------------------------------------------

!--- initialization of physics constants needed in the MYNN PBL scheme and already available from MPAS:
 cp      = con_cp
 cpv     = con_cpv
 cliq    = con_cliq
 cice    = con_cice
 ep_1    = con_ep1
 ep_2    = con_ep2
 grav    = con_grav
 karman  = con_karman
 p1000mb = con_p0
 r_d     = con_rd
 r_v     = con_rv
 rvovrd  = r_v/r_d
 svp1    = con_svp1
 svp2    = con_svp2
 svp3    = con_svp3
 svpt0   = con_svpt0
 xlf     = con_xlf
 xls     = con_xls
 xlv     = con_xlv

!--- initialization of derived physics constants needed in the MYNN PBL scheme:
 ep_3    = 1.-ep_2
 gtr     = grav/tref
 p608    = ep_1
 rcp     = r_d/cp
 t0c     = svpt0
 tv0     = p608*tref
 xlscp   = (xlv+xlf)/cp
 xlvcp   = xlv/cp

!ev      = xlv
!rk      = cp/r_d
!svp11   = svp1*1.e3
!tv1     = (1.+p608)*tref
!vk      = karman

 errmsg = " "
 errflg = 0

 end subroutine bl_mynn_init

!=================================================================================================================
 subroutine bl_mynn_finalize(errmsg,errflg)
!=================================================================================================================

 character(len=*),intent(out):: errmsg
 integer,intent(out):: errflg

!-----------------------------------------------------------------------------------------------------------------
!note: this subroutine currently does nothing.

 errmsg = ' '
 errflg = 0

 end subroutine bl_mynn_finalize

!=================================================================================================================
 subroutine bl_mynn_run &
            (initflag          , restart            , cycling            , &
             delt              , dz                 , dx                 , &
             znt               , u                  , v                  , &
             w                 , th                 , sqv                , &
             sqc               , sqi                , sqs                , &
             qnc               , qni                , qnwfa              , &
             qnifa             , qnbca              , qozone             , &
             p                 , exner              , rho                , &
             tt                , xland              , ts                 , &
             qsfc              , ps                 , ust                , &
             ch                , hfx                , qfx                , &
             rmol              , wspd               , uoce               , &
             voce              , qke                , qke_adv            , &
             tsq               , qsq                , cov                , &
             rublten           , rvblten            , rthblten           , &
             rqvblten          , rqcblten           , rqiblten           , &
             rqsblten          , rqncblten          , rqniblten          , &
             rqnwfablten       , rqnifablten        , rqnbcablten        , &
             rqozblten         , exch_h             , exch_m             , &
             pblh              , kpbl               , el_pbl             , &
             dqke              , qwt                , qshear             , &
             qbuoy             , qdiss              , sh                 , &
             sm                , qc_bl              , qi_bl              , &
             cldfra_bl         , icloud_bl          , bl_mynn_tkeadvect  , &
             bl_mynn_tkebudget , bl_mynn_cloudpdf   , bl_mynn_mixlength  , &
             bl_mynn_closure   , bl_mynn_stfunc     , bl_mynn_topdown    , &
             bl_mynn_edmf      , bl_mynn_edmf_dd    , bl_mynn_edmf_mom   , &
             bl_mynn_edmf_tke  , bl_mynn_mixscalars , bl_mynn_output     , &
             bl_mynn_cloudmix  , bl_mynn_mixqt      , bl_mynn_scaleaware , &
             bl_mynn_dheatopt  , edmf_a             , edmf_w             , &
             edmf_qt           , edmf_thl           , edmf_ent           , &
             edmf_qc           , sub_thl            , sub_sqv            , &
             det_thl           , det_sqv            , edmf_a_dd          , &
             edmf_w_dd         , edmf_qt_dd         , edmf_thl_dd        , &
             edmf_ent_dd       , edmf_qc_dd         , maxwidth           , &
             maxmf             , ztop_plume         , ktop_plume         , &
             spp_pbl           , pattern_spp_pbl    , rthraten           , &
             flag_qc           , flag_qi            , flag_qs            , &
             flag_qnc          , flag_qni           , flag_qnwfa         , &
             flag_qnifa        , flag_qnbca         , flag_qoz           , &
#if(WRF_CHEM == 1)
             mix_chem          , nchem              , kdvel              , &
             ndvel             , chem               , emis_ant_no        , &
             frp               , vdep               ,                      &
#endif
             its, ite , kts , kte , kme , errmsg , errflg                  &
            )    

!=================================================================================================================

!input arguments:
 logical,intent(in):: &
    flag_qc,flag_qi,flag_qs,flag_qoz,flag_qnc,flag_qni,flag_qnifa,flag_qnwfa,flag_qnbca

 logical,intent(in):: bl_mynn_edmf,bl_mynn_edmf_dd,bl_mynn_edmf_mom,bl_mynn_edmf_tke
 logical,intent(in):: bl_mynn_mixscalars,bl_mynn_cloudmix,bl_mynn_mixqt
 logical,intent(in):: bl_mynn_tkeadvect,bl_mynn_tkebudget
 logical,intent(in):: bl_mynn_output,bl_mynn_dheatopt,bl_mynn_scaleaware,bl_mynn_topdown

 logical,intent(in):: &
    restart,cycling

 integer,intent(in):: its,ite,kts,kte,kme

 integer,intent(in):: &
    initflag,icloud_bl,spp_pbl

 integer,intent(in):: &
    bl_mynn_cloudpdf,bl_mynn_mixlength,bl_mynn_stfunc
 
 real(kind=kind_phys),intent(in):: &
    bl_mynn_closure

 real(kind=kind_phys),intent(in):: &
    delt

 real(kind=kind_phys),intent(in),dimension(its:):: &
    dx,          &!
    xland,       &!
    ps,          &!
    ts,          &!
    qsfc,        &!
    ust,         &!
    ch,          &!
    hfx,         &!
    qfx,         &!
    rmol,        &!
    wspd,        &!
    uoce,        &!
    voce,        &!
    znt           !

 real(kind=kind_phys),intent(in),dimension(its:,:):: &
    dz,          &!
    u,           &!
    v,           &!
    th,          &!
    tt,          &!
    p,           &!
    exner,       &!
    rho,         &!
    rthraten      !

 real(kind=kind_phys),intent(in),dimension(its:,:):: &
    sqv,         &!
    sqc,         &!
    sqi,         &!
    sqs,         &!
    qnc,         &!
    qni,         &!
    qnifa,       &!
    qnwfa,       &!
    qnbca,       &!
    qozone        !

 real(kind=kind_phys),intent(in),dimension(its:,:):: &
    pattern_spp_pbl

 real(kind=kind_phys),intent(in),dimension(its:,:):: &
    w             !

!inout arguments:
 integer,intent(inout),dimension(its:):: &
    kpbl,        &!
    ktop_plume    !

 real(kind=kind_phys),intent(inout),dimension(its:):: &
    pblh

 real(kind=kind_phys),intent(inout),dimension(its:,:):: &
    cldfra_bl,   &!
    qc_bl,       &!
    qi_bl         !

 real(kind=kind_phys),intent(inout),dimension(its:,:):: &
    el_pbl,      &!
    qke,         &!
    qke_adv,     &!
    cov,         &!
    qsq,         &!
    tsq,         &!
    sh,          &!
    sm            !

 real(kind=kind_phys),intent(inout),dimension(its:,:):: &
    rublten,     &!
    rvblten,     &!
    rthblten,    &!
    rqvblten,    &!
    rqcblten,    &!
    rqiblten,    &!
    rqsblten,    &!
    rqncblten,   &!
    rqniblten,   &!
    rqnifablten, &!
    rqnwfablten, &!
    rqnbcablten, &!
    rqozblten     !

 real(kind=kind_phys),intent(inout),dimension(its:,:):: &
    edmf_a,      &!
    edmf_w,      &!
    edmf_qt,     &!
    edmf_thl,    &!
    edmf_ent,    &!
    edmf_qc,     &!
    sub_thl,     &!
    sub_sqv,     &!
    det_thl,     &!
    det_sqv       !

 real(kind=kind_phys),intent(inout),dimension(its:,:),optional:: &
    edmf_a_dd,   &!
    edmf_w_dd,   &!
    edmf_qt_dd,  &!
    edmf_thl_dd, &!
    edmf_ent_dd, &!
    edmf_qc_dd


!output arguments:
 character(len=*),intent(out):: &
    errmsg        ! output error message (-).

 integer,intent(out):: &
    errflg        ! output error flag (-).

 real(kind=kind_phys),intent(out),dimension(:):: &
    maxwidth,   &!
    maxmf,      &!
    ztop_plume

 real(kind=kind_phys),intent(out),dimension(its:,:):: &
    exch_h,      &!
    exch_m        !

 real(kind=kind_phys),intent(out),dimension(its:,:),optional:: &
    dqke,        &!
    qwt,         &!
    qshear,      &!
    qbuoy,       &!
    qdiss         !


!local variable and arrays:
 logical:: initialize_qke

 integer:: i,k

 real(kind=kind_phys):: qc_bl2,qi_bl2
 real(kind=kind_phys):: cpm,exnerg,flq,flqc,flqv,flt,fltv,phh,pmz,psig_bl,psig_shcu,sqcg,phi_m, &
                        th_sfc,zet,ts_decay

 real(kind=kind_phys),dimension(kts:kte):: cldfra_bl1_old,qc_bl1_old,qi_bl1_old
 real(kind=kind_phys),dimension(kts:kte):: qv1,qc1,qi1,qs1
 real(kind=kind_phys),dimension(kts:kte):: det_sqc,det_u,det_v,sub_u,sub_v

 real(kind=kind_phys),dimension(kts:kte):: pdc,pdk,pdq,pdt,sgm,sqw,thetav,thl,vq,vt,kzero

 real(kind=kind_phys),dimension(kts:kte):: dfh,dfm,dfq,qcd,tcd,diss_heat

 real(kind=kind_phys),dimension(kts:kte):: rstoch_col

 real(kind=kind_phys),dimension(kts:kte+1):: zw

 real(kind=kind_phys),dimension(kts:kte+1):: &
    s_aw1,s_awthl1,s_awqt1,s_awqv1,s_awqc1,s_awu1,s_awv1,s_awqke1,s_awqnc1,s_awqni1, &
    s_awqnwfa1,s_awqnifa1,s_awqnbca1

 real(kind=kind_phys),dimension(kts:kte+1):: &
    sd_aw1,sd_awthl1,sd_awqt1,sd_awqv1,sd_awqc1,sd_awu1,sd_awv1,sd_awqke1

!JOE-top-down diffusion
 logical :: cloudflg
 integer :: kk,kminrad

 real(kind=kind_phys),parameter:: pfac =2.0, zfmin = 0.01, phifac=8.0
 real(kind=kind_phys):: maxkhtopdown
 real(kind=kind_phys):: bfxpbl,dthvx,tmp1,temps,templ,zl1,wstar3_2
 real(kind=kind_phys):: ent_eff,radsum,radflux,we,rcldb,rvls,minrad,zminrad
 real(kind=kind_phys),dimension(kts:kte):: khtopdown,zfac,wscalek2,zfacent,tkeprodtd
!JOE-end top down


!local 1D input arguments:
 real(kind=kind_phys):: dx1,xland1,ps1,ts1,qsfc1,ust1,ch1,hfx1,qfx1,rmol1,wspd1, &
    uoce1,voce1,znt1
 real(kind=kind_phys),dimension(kts:kte):: &
    dz1,u1,v1,th1,tk1,p1,ex1,rho1,qnc1,qni1,qnifa1,qnwfa1,qnbca1,qozone1,rthraten1,sqv1,sqc1,sqi1,sqs1
 real(kind=kind_phys),dimension(kts:kme):: w1

!local 1D inout arguments:
 integer:: kpbl1,ktop_plume1

 real(kind=kind_phys):: pblh1
 real(kind=kind_phys),dimension(kts:kte):: cldfra_bl1,qc_bl1,qi_bl1
 real(kind=kind_phys),dimension(kts:kte):: el_pbl1,qke1,qke_adv1,cov1,qsq1,tsq1,sh1,sm1
 real(kind=kind_phys),dimension(kts:kte):: du1,dv1,dth1,dqv1,dqc1,dqi1,dqs1,dqnc1,dqni1,dqnifa1,dqnwfa1, &
                                           dqnbca1,dqozone1
 real(kind=kind_phys),dimension(kts:kte):: edmf_a1,edmf_w1,edmf_qt1,edmf_thl1,edmf_ent1,edmf_qc1,sub_thl1, &
                                           sub_sqv1,det_thl1,det_sqv1
 real(kind=kind_phys),dimension(kts:kte):: edmf_a_dd1,edmf_w_dd1,edmf_qt_dd1,edmf_thl_dd1,edmf_ent_dd1, &
                                           edmf_qc_dd1

!local 1D output arguments:
 real(kind=kind_phys):: maxwidth1,maxmf1,ztop_plume1
 real(kind=kind_phys),dimension(kts:kte):: exch_h1,exch_m1

 real(kind=kind_phys),dimension(kts:kte):: dqke1,qwt1,qshear1,qbuoy1,qdiss1

!substepping TKE:
 integer:: nsub
 real(kind=kind_phys):: delt2



!VARIABLES NEEDED FOR MIXING OF CHEMICAL SPECIES:
#if(WRF_CHEM == 1)
!--- inputs:
 logical,intent(in):: mix_chem
 integer,intent(in):: nchem,kdvel,ndvel
 real(kind=kind_phys),intent(in),dimension(its:ite),optional:: frp,emis_ant_no
 real(kind=kind_phys),intent(in),dimension(its:ite,ndvel):: vdep

!--- inouts:
 real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte,nchem):: chem
#else
 logical,parameter:: mix_chem = .false.
 integer,parameter:: nchem = 1
 integer,parameter:: kdvel = 1
 integer,parameter:: ndvel = 1
#endif
!--- local variables and arrays:
 logical,parameter:: rrfs_sd   = .false.
 logical,parameter:: smoke_dbg = .false.
 logical,parameter:: enh_mix   = .false.

 integer:: ic 
 real(kind=kind_phys):: emis_ant_no1,frp1
 real(kind=kind_phys),dimension(ndvel):: vd1
 real(kind=kind_phys),dimension(kts:kte,nchem)::   chem1
 real(kind=kind_phys),dimension(kts:kte+1,nchem):: s_awchem1          
 !END VARIABLES NEEDED FOR MIXING OF CHEMICAL SPECIES.

!-----------------------------------------------------------------------------------------------------------------

 errmsg = " "
 errflg = 0

 do i = its,ite

    if(present(dqke)) then
       do k = kts,kte
          dqke(i,k) = qke(i,k)
       enddo
    endif

!--- initialization of 2D inout tendencies:
    do k = kts,kte
       rublten(i,k)     = 0._kind_phys
       rvblten(i,k)     = 0._kind_phys
       rthblten(i,k)    = 0._kind_phys
       rqvblten(i,k)    = 0._kind_phys
       rqcblten(i,k)    = 0._kind_phys       
       rqiblten(i,k)    = 0._kind_phys       
       rqsblten(i,k)    = 0._kind_phys
       rqncblten(i,k)   = 0._kind_phys       
       rqniblten(i,k)   = 0._kind_phys       
       rqnifablten(i,k) = 0._kind_phys       
       rqnwfablten(i,k) = 0._kind_phys       
       rqnbcablten(i,k) = 0._kind_phys
       rqozblten(i,k)   = 0._kind_phys
    enddo

!--- initialization of 2D output variables:
    ktop_plume(i) = 0
    maxwidth(i)   = 0._kind_phys
    maxmf(i)      = 0._kind_phys
    ztop_plume(i) = 0._kind_phys

!--- initialization of 1D input variables using 2D input variables:
    dx1    = dx(i)
    xland1 = xland(i)
    ps1    = ps(i)
    ts1    = ts(i)
    qsfc1  = qsfc(i)
    ust1   = ust(i)
    ch1    = ch(i)
    hfx1   = hfx(i)
    qfx1   = qfx(i)
    rmol1  = rmol(i)
    wspd1  = wspd(i)
    uoce1  = uoce(i)
    voce1  = voce(i)
    znt1   = znt(i)

    do k = kts,kte
       dz1(k)        = dz(i,k)
       u1(k)         = u(i,k)
       v1(k)         = v(i,k)
       w1(k)         = w(i,k)
       th1(k)        = th(i,k)
       tk1(k)        = tt(i,k)
       p1(k)         = p(i,k)
       ex1(k)        = exner(i,k)
       rho1(k)       = rho(i,k)
       sh1(k)        = sh(i,k)
       sm1(k)        = sm(i,k)
       rthraten1(k)  = rthraten(i,k)
       sqv1(k)       = sqv(i,k)
       sqc1(k)       = sqc(i,k)
       sqi1(k)       = sqi(i,k)
       sqs1(k)       = sqs(i,k)
       qnc1(k)       = qnc(i,k)
       qni1(k)       = qni(i,k)
       qnifa1(k)     = qnifa(i,k)
       qnwfa1(k)     = qnwfa(i,k)
       qnbca1(k)     = qnbca(i,k)
       qozone1(k)    = qozone(i,k)
       kzero(k)      = 0._kind_phys
    enddo
    do k = kte,kte+1
       w1(k) = w(i,k)
    enddo

!--- initialization of the PBL stochastic forcing:
    if(spp_pbl .eq. 1) then
       do k = kts,kte
          rstoch_col(k) = pattern_spp_pbl(i,k)
       enddo
    else
       do k = kts,kte
          rstoch_col(k) = 0._kind_phys
       enddo
    endif


!--- initialization of 1D inout variables using 2D inout variables:
    kpbl1 = kpbl(i)
    pblh1 = pblh(i)

    do k = kts,kte
       cldfra_bl1(k) = cldfra_bl(i,k)
       qc_bl1(k)     = qc_bl(i,k)
       qi_bl1(k)     = qi_bl(i,k)
    enddo

    do k = kts,kte
       el_pbl1(k)  = el_pbl(i,k)
       qke1(k)     = qke(i,k)
       qke_adv1(k) = qke_adv(i,k)
       cov1(k)     = cov(i,k)
       qsq1(k)     = qsq(i,k)
       tsq1(k)     = tsq(i,k)
       sh1(k)      = sh(i,k)
       sm1(k)      = sm(i,k)
    enddo

!--- initialization of 1D local variables:
    ktop_plume1  = 0
    maxwidth1    = 0._kind_phys
    maxmf1       = 0._kind_phys
    ztop_plume1  = 0._kind_phys
    maxkhtopdown = 0._kind_phys

    do k = kts,kte
       du1(k)      = 0._kind_phys
       dv1(k)      = 0._kind_phys
       dth1(k)     = 0._kind_phys
       dqv1(k)     = 0._kind_phys
       dqc1(k)     = 0._kind_phys
       dqi1(k)     = 0._kind_phys
       dqs1(k)     = 0._kind_phys
       dqnc1(k)    = 0._kind_phys
       dqni1(k)    = 0._kind_phys
       dqnifa1(k)  = 0._kind_phys
       dqnwfa1(k)  = 0._kind_phys
       dqnbca1(k)  = 0._kind_phys
       dqozone1(k) = 0._kind_phys
    enddo
    do k = kts,kte
       edmf_a1(k)   = 0._kind_phys
       edmf_w1(k)   = 0._kind_phys
       edmf_qc1(k)  = 0._kind_phys
       edmf_ent1(k) = 0._kind_phys
       edmf_qt1(k)  = 0._kind_phys
       edmf_thl1(k) = 0._kind_phys
       sub_thl1(k)  = 0._kind_phys
       sub_sqv1(k)  = 0._kind_phys
       det_thl1(k)  = 0._kind_phys
       det_sqv1(k)  = 0._kind_phys

       edmf_a_dd1(k)   = 0._kind_phys
       edmf_w_dd1(k)   = 0._kind_phys
       edmf_qc_dd1(k)  = 0._kind_phys
       edmf_ent_dd1(k) = 0._kind_phys
       edmf_qt_dd1(k)  = 0._kind_phys
       edmf_thl_dd1(k) = 0._kind_phys
    enddo
    do k = kts,kte
       dqke1(k)   = 0._kind_phys
       qwt1(k)    = 0._kind_phys
       qshear1(k) = 0._kind_phys
       qbuoy1(k)  = 0._kind_phys
       qdiss1(k)  = 0._kind_phys
       exch_h1(k) = 0._kind_phys
       exch_m1(k) = 0._kind_phys
    enddo
    do k = kts,kte
       sub_u(k)   = 0._kind_phys
       sub_v(k)   = 0._kind_phys
       det_sqc(k) = 0._kind_phys
       det_u(k)   = 0._kind_phys
       det_v(k)   = 0._kind_phys
    enddo
    do k = kts,kte+1
       s_aw1(k)      = 0._kind_phys
       s_awthl1(k)   = 0._kind_phys
       s_awqt1(k)    = 0._kind_phys
       s_awqv1(k)    = 0._kind_phys
       s_awqc1(k)    = 0._kind_phys
       s_awu1(k)     = 0._kind_phys
       s_awv1(k)     = 0._kind_phys
       s_awqke1(k)   = 0._kind_phys
       s_awqnc1(k)   = 0._kind_phys
       s_awqni1(k)   = 0._kind_phys
       s_awqnwfa1(k) = 0._kind_phys
       s_awqnifa1(k) = 0._kind_phys
       s_awqnbca1(k) = 0._kind_phys
    enddo
    do k = kts,kte+1
       sd_aw1(k)     = 0._kind_phys
       sd_awthl1(k)  = 0._kind_phys
       sd_awqt1(k)   = 0._kind_phys
       sd_awqv1(k)   = 0._kind_phys
       sd_awqc1(k)   = 0._kind_phys
       sd_awu1(k)    = 0._kind_phys
       sd_awv1(k)    = 0._kind_phys
       sd_awqke1(k)  = 0._kind_phys
    enddo
    do k = kts,kte
       cldfra_bl1_old(k) = 0._kind_phys
       qc_bl1_old(k)     = 0._kind_phys
       qi_bl1_old(k)     = 0._kind_phys
    enddo

    do k = kts,kte
       qv1(k) = sqv1(k)/(1.-sqv1(k))
       qc1(k) = sqc1(k)/(1.-sqv1(k))
       qi1(k) = sqi1(k)/(1.-sqv1(k))
       qs1(k) = sqs1(k)/(1.-sqv1(k))
    enddo

    k = kts
    zw(k) = 0._kind_phys
    do k = kts+1,kte+1
       zw(k) = zw(k-1) + dz1(k-1)
    enddo

!INITIALIZATION OF LOCAL CHEMICAL SPECIES:
#if(WRF_CHEM == 1)
    do ic = 1,nchem
       vd1(ic) = vdep(i,ic)
       do k = kts,kte
          chem1(k,ic) = chem(i,k,ic)
       enddo
    enddo
    if(present(emis_ant_no) .and. present(frp)) then
       emis_ant_no1 = emis_ant_no(i)
       frp1 = frp(i)
    else
       emis_ant_no1 = 0._kind_phys
       frp1 = 0._kind_phys
    endif
    !END INITIALIZATION OF LOCAL CHEMICAL SPECIES.
#else
    do ic = 1,nchem
       vd1(ic) = 0._kind_phys
       do k = kts,kte
          chem1(k,ic) = 0._kind_phys
       enddo
    enddo
    emis_ant_no1 = 0._kind_phys
    frp1 = 0._kind_phys
#endif
    do ic = 1,nchem
       do k = kts,kte+1
          s_awchem1(k,ic) = 0._kind_phys
       enddo
    enddo
!END INITIALIZATION OF LOCAL CHEMICAL SPECIES.


    do k = kts,kte
       !keep snow out for now - increase ceiling bias
       sqw(k) = sqv1(k)+sqc1(k)+sqi1(k) !+sqs1(k)
       thl(k) = th1(k) - xlvcp/ex1(k)*sqc1(k) - xlscp/ex1(k)*(sqi1(k))!+sqs1(k))
       thetav(k) = th1(k)*(1.+0.608*sqv1(k))

       !Use form from Tripoli and Cotton (1981) with their
       !suggested min temperature to improve accuracy.    
       !thl(k)=th(k)*(1.- xlvcp/MAX(tk1(k),TKmin)*sqc1(k) - xlscp/MAX(tk1(k),TKmin)*sqi1(k))
       !thetav(k) = th1(k)*(1.+p608)*sqv1(k)
    enddo

!-----------------------------------------------------------------------------------------------------------------
!initflag > 0:
!-----------------------------------------------------------------------------------------------------------------
    if(initflag > 0 .and. .not.restart) then

       !test to see if we want to initialize qke1:
       if((restart .or. cycling)) then
          if(qke1(kts) < 0.0002) then
             initialize_qke = .true.
          else
             initialize_qke = .false.
          endif
       else ! not cycling or restarting:
          initialize_qke = .true.
       endif

       if(.not.restart .or. .not.cycling) then
          do k = kts,kte
             sh1(k)        = 0._kind_phys
             sm1(k)        = 0._kind_phys
             el_pbl1(k)    = 0._kind_phys
             tsq1(k)       = 0._kind_phys
             qsq1(k)       = 0._kind_phys
             cov1(k)       = 0._kind_phys
             cldfra_bl1(k) = 0._kind_phys
             qc_bl1(k)     = 0._kind_phys
             qi_bl1(k)     = 0._kind_phys
             qke1(k)       = 0._kind_phys
          enddo
       endif
       do k = kts,kte
          cldfra_bl1_old(k) = 0._kind_phys
          qc_bl1_old(k)     = 0._kind_phys
          qi_bl1_old(k)     = 0._kind_phys
       enddo

       if(initialize_qke) then
          do k = kts,kte
             qke1(k)=5.*ust1*max((ust1*700.-zw(k))/(max(ust1,0.01)*700.),0.01)
          enddo
       endif

       !--- computes the PBL height:
       call get_pblh(kts,kte,pblh1,thetav,qke1,zw,dz1,xland1,kpbl1)

       !--- computes the similarity functions:
       if(bl_mynn_scaleaware) then
          call scale_aware(dx1,pblh1,psig_bl,psig_shcu)
       else
          psig_bl = 1._kind_phys
          psig_shcu = 1._kind_phys
       endif

       !--- calls mym_initialize:
       call mym_initialize(                   & 
                kts,kte,xland1,               &
                dz1,dx1,zw,                   &
                u1,v1,thl,sqv1,               &
                pblh1,th1,thetav,sh1,sm1,     &
                ust1, rmol1,                  &
                el_pbl1,qke1,tsq1,qsq1,cov1,  &
                psig_bl,cldfra_bl1,           &
                bl_mynn_mixlength,            &
                edmf_w1,edmf_a1,              &
                initialize_qke,               &
                spp_pbl,rstoch_col)

    endif
!-----------------------------------------------------------------------------------------------------------------
!end initflag > 0:
!-----------------------------------------------------------------------------------------------------------------


    if(bl_mynn_tkeadvect) then
       do k = kts,kte
          qke1(k) = qke_adv1(k)
       enddo
    endif
    !Joe-TKE budget:
    if(bl_mynn_tkebudget) then
       do k = kts,kte
          dqke1(k) = qke1(k)
       enddo
    endif
    if(icloud_bl > 0) then
       do k = kts,kte
          cldfra_bl1_old(k) = cldfra_bl1(k)
          qc_bl1_old(k)     = qc_bl1(k)
          qi_bl1_old(k)     = qi_bl1(k)
       enddo
    endif

    !--- computes the PBL height:
    call get_pblh(kts,kte,pblh1,thetav,qke1,zw,dz1,xland1,kpbl1)

    !--- computes the similarity functions:
    if(bl_mynn_scaleaware) then
       call scale_aware(dx1,pblh1,psig_bl,psig_shcu)
    else
       psig_bl   = 1._kind_phys
       psig_shcu = 1._kind_phys
    endif

    sqcg   = 0.0   !ill-defined variable; qcg has been removed
    cpm    = cp*(1.+0.84*qv1(kts))
    exnerg = (ps1/p1000mb)**rcp

    !-----------------------------------------------------
    !ORIGINAL CODE
    !flt = hfx(i)/( rho(i,kts)*cpm ) &
    ! +xlvcp*ch(i)*(sqc(kts)/exner(i,kts) -sqcg/exnerg)
    !flq = qfx(i)/  rho(i,kts)       &
    !    -ch(i)*(sqc(kts)   -sqcg )
    !-----------------------------------------------------
    flqv   = qfx1/rho1(kts)
    flqc   = 0.0 !currently no sea-spray fluxes, fog settling hangled elsewhere
    th_sfc = ts1/ex1(kts)

    !--- turbulent flux for the TKE voundary conditions:
    flq  = flqv + flqc                                  ! Latent
    flt  = hfx1/(rho1(kts)*cpm ) - xlvcp*flqc/ex1(kts)  ! Temperature flux
    fltv = flt + flqv*p608*th_sfc                       ! Virtual temperature flux  

    !--- update 1/L using updated sfc heat flux and friction velocity:
    rmol1 = -karman*gtr*fltv/max(ust1**3,1.0e-6)
    zet = 0.5*dz1(kts)*rmol1
    zet = max(zet, -20.)
    zet = min(zet,  20.)                                         

    !if(i.eq.idbg)print*,"updated z/L=",zet                      
    if(bl_mynn_stfunc == 0) then                                
       !original Kansas-type stability functions:
       if(zet >= 0.0) then                                    
          pmz = 1.0 + (cphm_st-1.0) * zet
          phh = 1.0 +  cphh_st      * zet
       else
          pmz = 1.0/    (1.0-cphm_unst*zet)**0.25 - zet          
          phh = 1.0/sqrt(1.0-cphh_unst*zet)                      
       endif
       phi_m = pmz + zet
    else
       !updated stability functions (Puhales, 2020):
       phi_m = phim(zet)
       pmz   = phi_m - zet                                       
       phh   = phih(zet)                                         
    endif 

    !call mym_condensation() to calculate the nonconvective component of the subgrid-scale cloud fraction
    !and mixing ratio as well as the functions used to calculate the buoyancy flux. Different cloud PDFs
    !can be selected by use of the namelist parameter bl_mynn_cloudpdf:
    do k = kts,kte
       vt(k)  = 0._kind_phys
       vq(k)  = 0._kind_phys
       sgm(k) = 0._kind_phys
    enddo

    call  mym_condensation(kts,kte,          &
              dx1,dz1,zw,xland1,             &
              thl,sqw,sqv1,sqc1,sqi1,sqs1,   &
              p1,ex1,tsq1,qsq1,cov1,         &
              sh1,el_pbl1,bl_mynn_cloudpdf,  &
              qc_bl1,qi_bl1,cldfra_bl1,      &
              pblh1,hfx1,                    &
              vt,vq,th1,sgm,rmol1,           &
              spp_pbl,rstoch_col)


    !add TKE source driven by cloud top cooling. calculate the buoyancy production of tke from cloud-top
    !cooling when bl_mynn_topdown = .true.
    if(bl_mynn_topdown)then
       call topdown_cloudrad(kts,kte,dz1,zw,fltv,             &
                    xland1,kpbl1,pblh1,                       &
                    sqc1,sqi1,sqw,thl,th1,ex1,p1,rho1,thetav, &
                    cldfra_bl1,rthraten1,                     &
                    maxkhtopdown,khtopdown,tkeprodtd)
    else
       maxkhtopdown = 0._kind_phys
       do k = kts,kte
          khtopdown(k) = 0._kind_phys
          tkeprodtd(k) = 0._kind_phys
       enddo
    endif


    !--- calls subroutine dmp_mf():
    if(bl_mynn_edmf) then
       call dmp_mf( i,                             &
                kts,kte,delt,zw,dz1,p1,rho1,       &
                bl_mynn_edmf_mom,                  &
                bl_mynn_edmf_tke,                  &
                bl_mynn_mixscalars,                &
                u1,v1,w1,th1,thl,thetav,tk1,       &
                sqw,sqv1,sqc1,qke1,                &
                qnc1,qni1,qnwfa1,qnifa1,qnbca1,    &
                ex1,vt,vq,sgm,                     &
                ust1,flt,fltv,flq,flqv,            &
                pblh1,kpbl1,dx1,                   &
                xland1,th_sfc,                     &
                !now outputs - tendencies
                !dth1mf,dqv1mf,dqc1mf,du1mf,dv1mf, &
                !outputs - updraft properties
                edmf_a1,edmf_w1,edmf_qt1,          &
                edmf_thl1,edmf_ent1,edmf_qc1,      &
                !for the solver
                s_aw1,s_awthl1,s_awqt1,            &
                s_awqv1,s_awqc1,                   &
                s_awu1,s_awv1,s_awqke1,            &
                s_awqnc1,s_awqni1,                 &
                s_awqnwfa1,s_awqnifa1,s_awqnbca1,  &
                sub_thl1,sub_sqv1,                 &
                sub_u,sub_v,                       &
                det_thl1,det_sqv1,det_sqc,         &
                det_u,det_v,                       &
                !chem/smoke mixing
                nchem,chem1,s_awchem1,             &
                mix_chem,                          &
                qc_bl1,cldfra_bl1,                 &
                qc_bl1_old,cldfra_bl1_old,         &
                flag_qc,flag_qi,                   &
                flag_qnc,flag_qni,                 &
                flag_qnwfa,flag_qnifa,flag_qnbca,  &
                psig_shcu,                         &
                maxwidth1,ktop_plume1,             &
                maxmf1,ztop_plume1,                &
                spp_pbl,rstoch_col)

       if(bl_mynn_edmf_dd) then
          call ddmf_jpl(kts,kte,delt,zw,dz1,p1,         &
                    u1,v1,th1,thl,thetav,tk1,           &
                    sqw,sqv1,sqc1,rho1,ex1,             &
                    ust1,flt,flq,                       &
                    pblh1,kpbl1,                        &
                    edmf_a_dd1,edmf_w_dd1,edmf_qt_dd1,  &
                    edmf_thl_dd1,edmf_ent_dd1,          &
                    edmf_qc_dd1,                        &
                    sd_aw1,sd_awthl1,sd_awqt1,          &
                    sd_awqv1,sd_awqc1,sd_awu1,sd_awv1,  &
                    sd_awqke1,                          &
                    qc_bl1,cldfra_bl1,                  &
                    rthraten)
       endif
    endif


    !--- capability to substep the eddy-diffusivity portion:
    !do nsub = 1,2
    delt2 = delt !*0.5    !only works if topdown=0

    call mym_turbulence                        &
             (kts,kte,xland1,bl_mynn_closure,  &
             dz1,dx1,zw,                       &
             u1,v1,thl,thetav,sqc1,sqw,        &
             qke1,tsq1,qsq1,cov1,              &
             vt,vq,                            &
             rmol1,flt,fltv,flq,               &
             pblh1,th1,                        &
             sh1,sm1,el_pbl1,                  &
             dfm,dfh,dfq,                      &
             tcd,qcd,pdk,                      &
             pdt,pdq,pdc,                      &
             qwt1,qshear1,qbuoy1,qdiss1,       &
             bl_mynn_tkebudget,                &
             psig_bl,psig_shcu,                &
             cldfra_bl1,bl_mynn_mixlength,     &
             edmf_w1,edmf_a1,                  &
             tkeprodtd,                        &
             spp_pbl,rstoch_col)


    !--- calls subroutine mym_predict() to solve TKE:
    call mym_predict                             &
             (kts,kte,bl_mynn_closure,           &
              delt2,dz1,                         &
              ust1,flt,flq,pmz,phh,              &
              el_pbl1,dfq,rho1,pdk,pdt,pdq,pdc,  &
              qke1,tsq1,qsq1,cov1,               &
              s_aw1,s_awqke1,bl_mynn_edmf_tke,   &
              qwt1,qdiss1,bl_mynn_tkebudget)      ! TKE budget (Puhales 2020)

    if(bl_mynn_dheatopt) then
       do k = kts,kte-1
          !set max dissipative heating rate to 7.2 K per hour
          diss_heat(k) = min(max(1.0*(qke1(k)**1.5)/(b1*max(0.5*(el_pbl1(k)+el_pbl1(k+1)),1.))/cp,0.0),0.002)

          !limit heating above 100 mb:
          diss_heat(k) = diss_heat(k) * exp(-10000./MAX(p1(k),1.))
       enddo
       diss_heat(kte) = 0.
    else
       do k = kts,kte
          diss_heat(k) = 0.
       enddo
    endif       


    !--- call to subroutine mynn_tendencies:
    call mynn_tendencies(kts,kte, &
              delt,dz1,rho1,                     &
              u1,v1,th1,tk1,qv1,                 &
              qc1,qi1,kzero,qnc1,qni1,           & !kzero replaces qs1 - not mixing snow
              ps1,p1,ex1,thl,                    &
              sqv1,sqc1,sqi1,kzero,sqw,          & !kzero replaces sqs - not mxing snow
              qnwfa1,qnifa1,qnbca1,qozone1,      &
              ust1,flt,flq,flqv,flqc,            &
              wspd1,uoce1,voce1,                 &
              tsq1,qsq1,cov1,                    &
              tcd,qcd,                           &
              dfm,dfh,dfq,                       &
              du1,dv1,dth1,dqv1,                 &
              dqc1,dqi1,dqs1,dqnc1,dqni1,        &
              dqnwfa1,dqnifa1,dqnbca1,           &
              dqozone1,                          &
              diss_heat,                         &
              !mass flux components
              s_aw1,s_awthl1,s_awqt1,            &
              s_awqv1,s_awqc1,s_awu1,s_awv1,     &
              s_awqnc1,s_awqni1,                 &
              s_awqnwfa1,s_awqnifa1,s_awqnbca1,  &
              sd_aw1,sd_awthl1,sd_awqt1,         &
              sd_awqv1,sd_awqc1,                 &
              sd_awu1,sd_awv1,                   &
              sub_thl1,sub_sqv1,                 &
              sub_u,sub_v,                       &
              det_thl1,det_sqv1,det_sqc,         &
              det_u,det_v,                       &
              flag_qc,flag_qi,flag_qnc,          &
              flag_qni,flag_qs,                  &
              flag_qnwfa,flag_qnifa,             &
              flag_qnbca,flag_qoz,               &
              cldfra_bl1,                        &
              bl_mynn_cloudmix,                  &
              bl_mynn_mixqt,                     &
              bl_mynn_edmf_mom,                  &
              bl_mynn_mixscalars)


    !--- call to subroutine mynn_mix_chem for PBL and tropospheric mixing of
    !    chemical species:
    if(mix_chem) then
       if(rrfs_sd) then
          call mynn_mix_chem(kts,kte,   &
                    delt,dz1,pblh1,     &
                    nchem,kdvel,ndvel,  &
                    chem1,vd1,          &
                    rho1,flt,           &
                    tcd,qcd,            &
                    dfh,                &
                    s_aw1,s_awchem1,    &
                    emis_ant_no1,       &
                    frp1,rrfs_sd,       &
                    enh_mix,smoke_dbg)
       else
          call mynn_mix_chem(kts,kte,   &
                    delt,dz1,pblh1,     &
                    nchem,kdvel,ndvel,  &
                    chem1,vd1,          &
                    rho1,flt,           &
                    tcd,qcd,            &
                    dfh,                &
                    s_aw1,s_awchem1,    &
                    zero,               &
                    zero,rrfs_sd,       &
                    enh_mix,smoke_dbg)
       endif
    endif
#if(WRF == 1)
    !directly updates chem3 instead of computing a tendency:
    do ic = 1,nchem
       do k = kts,kte
          chem(i,k,ic) = max(1.e-12,chem1(k,ic))
       enddo
    enddo
#endif


    !--- computes the exchange coefficients:
    call retrieve_exchange_coeffs(kts,kte,dfm,dfh,dz1,exch_m1,exch_h1)


!-----------------------------------------------------------------------------------------------------------------
!begin output of 2D variables:
!-----------------------------------------------------------------------------------------------------------------
    !output tendencies:
    do k = kts,kte
       rublten(i,k)  = du1(k)
       rvblten(i,k)  = dv1(k)
       rthblten(i,k) = dth1(k)
       rqvblten(i,k) = dqv1(k)
    enddo
    if(bl_mynn_cloudmix .and. flag_qc) then       
       do k = kts,kte
          rqcblten(i,k) = dqc1(k)
       enddo
    endif
    if(bl_mynn_cloudmix .and. flag_qi) then
       do k = kts,kte
         rqiblten(i,k) = dqi1(k)
       enddo
    endif
    if(bl_mynn_cloudmix .and. flag_qs) then
       do k = kts,kte
         rqsblten(i,k) = dqs1(k)
       enddo
    endif
    if(bl_mynn_cloudmix .and. bl_mynn_mixscalars .and. flag_qnc) then
       do k = kts,kte
          rqncblten(i,k) = dqnc1(k)
       enddo
    endif
    if(bl_mynn_cloudmix .and. bl_mynn_mixscalars .and. flag_qni) then
       do k = kts,kte
          rqniblten(i,k) = dqni1(k)
       enddo
    endif
    if(bl_mynn_cloudmix .and. bl_mynn_mixscalars .and. flag_qnifa) then
       do k = kts,kte
          rqnifablten(i,k) = dqnifa1(k)
       enddo
    endif
    if(bl_mynn_cloudmix .and. bl_mynn_mixscalars .and. flag_qnwfa) then
       do k = kts,kte
          rqnwfablten(i,k) = dqnwfa1(k)
       enddo
    endif
    if(bl_mynn_cloudmix .and. bl_mynn_mixscalars .and. flag_qnbca) then
       do k = kts,kte
          rqnbcablten(i,k) = dqnbca1(k)
       enddo
    endif
    do k = kts,kte
       rqozblten(i,k) = 0._kind_phys
    enddo

    !inout arrays:
    kpbl(i)       = kpbl1
    ktop_plume(i) = ktop_plume1

    pblh(i)  = pblh1

    do k = kts,kte
       cldfra_bl(i,k) = cldfra_bl1(k)
       qc_bl(i,k)     = qc_bl1(k)
       qi_bl(i,k)     = qi_bl1(k)
    enddo

    do k = kts,kte
       el_pbl(i,k)  = el_pbl1(k)
       qke(i,k)     = qke1(k)
       qke_adv(i,k) = qke_adv1(k)
       cov(i,k)     = cov1(k)
       qsq(i,k)     = qsq1(k)
       tsq(i,k)     = tsq1(k)
       sh(i,k)      = sh1(k)
       sm(i,k)      = sm1(k)
    enddo


    !the TKE budget is now given in m**2/s**-3 (Puhales, 2020):
    if(present(qwt)   .and. present(qbuoy) .and. present(qshear) .and. &
       present(qdiss) .and. present(dqke)) then
       if(bl_mynn_tkebudget) then
          !lower boundary condtions (using similarity relationships such as the prognostic equation for Qke)
          k = kts
          qshear1(k) = 4.*(ust1**3*phi_m/(karman*dz1(k)))-qshear1(k+1) ! staggered
          qbuoy1(k)  = 4.*(-ust1**3*zet/(karman*dz1(k)))-qbuoy1(k+1)   ! staggered

          !unstaggering shear and buoy and trasnfering all TKE budget to 2D arrays:
          do k = kts,kte-1
             qshear(i,k) = 0.5*(qshear1(k) + qshear1(k+1)) ! unstaggering in z
             qbuoy(i,k)  = 0.5*(qbuoy1(k)  + qbuoy1(k+1))  ! unstaggering in z
             qwt(i,k)    = qwt1(k)
             qdiss(i,k)  = qdiss1(k)
             dqke(i,k)   = (qke1(k)-dqke(i,k))*0.5/delt
          enddo
          !upper boundary conditions
          k = kte
          qshear(i,k) = 0._kind_phys
          qbuoy(i,k)  = 0._kind_phys
          qwt(i,k)    = 0._kind_phys
          qdiss(i,k)  = 0._kind_phys
          dqke(i,k)   = 0._kind_phys
       else
          do k = kts,kte
             qshear(i,k) = 0._kind_phys
             qbuoy(i,k)  = 0._kind_phys
             qwt(i,k)    = 0._kind_phys
             qdiss(i,k)  = 0._kind_phys
             dqke(i,k)   = 0._kind_phys
          enddo
       endif
    endif


    !optional inout arrays for updraft/downdraft properties:
    if(bl_mynn_edmf .and. bl_mynn_output) then
       do k = kts,kte
          edmf_a(i,k)   = edmf_a1(k)
          edmf_w(i,k)   = edmf_w1(k)
          edmf_qt(i,k)  = edmf_qt1(k)
          edmf_thl(i,k) = edmf_thl1(k)
          edmf_ent(i,k) = edmf_ent1(k)
          edmf_qc(i,k)  = edmf_qc1(k)
          sub_thl(i,k)  = sub_thl1(k)
          sub_sqv(i,k)  = sub_sqv1(k)
          det_thl(i,k)  = det_thl1(k)
          det_sqv(i,k)  = det_sqv1(k)
       enddo
    else
       do k = kts,kte
          edmf_a(i,k)   = 0._kind_phys
          edmf_w(i,k)   = 0._kind_phys
          edmf_qt(i,k)  = 0._kind_phys
          edmf_thl(i,k) = 0._kind_phys
          edmf_ent(i,k) = 0._kind_phys
          edmf_qc(i,k)  = 0._kind_phys
          sub_thl(i,k)  = 0._kind_phys
          sub_sqv(i,k)  = 0._kind_phys
          det_thl(i,k)  = 0._kind_phys
          det_sqv(i,k)  = 0._kind_phys
       enddo
    endif
    if(bl_mynn_edmf_dd .and. bl_mynn_output) then
       if(present(edmf_a_dd)   .and. present(edmf_w_dd)   .and. present(edmf_qt_dd) .and. &
          present(edmf_thl_dd) .and. present(edmf_ent_dd) .and. present(edmf_qc_dd)) then
          do k = kts,kte
             edmf_a_dd(i,k)   = edmf_a_dd1(k)
             edmf_w_dd(i,k)   = edmf_w_dd1(k)
             edmf_qt_dd(i,k)  = edmf_qt_dd1(k)
             edmf_thl_dd(i,k) = edmf_thl_dd1(k)
             edmf_ent_dd(i,k) = edmf_ent_dd1(k)
             edmf_qc_dd(i,k)  = edmf_qc_dd1(k)
          enddo
       endif
    endif

    !output arrays:
    maxwidth(i)   = maxwidth1
    maxmf(i)      = maxmf1
    ztop_plume(i) = ztop_plume1

    do k = kts,kte
       exch_h(i,k) = exch_h1(k)
       exch_m(i,k) = exch_m1(k)
    enddo

 enddo

 end subroutine bl_mynn_run

!=================================================================================================================
 end module bl_mynn
!=================================================================================================================
